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VESSEL SEGMENTATION USING VESSELNESS AND EDGENESS 
CROSS-REFERENCE TO RELATED APPLICATION 

This application claims priority to U.S*. Provisional Patent Application Serial 
No. 60/525,603, filed November 26, 2003, which is incorporated herein by reference 
in its entirety. 

BACKGROUND 

The present disclosure relates to the performance of virtual examinations. 
More particularly, the disclosure provides a system and method for vessel 
segmentation. 

Two-dimensional (2D) visualization of human organs using medical imaging 
devices has been widely used for patient diagnosis. Currently available medical 
imaging devices include computed tomography (CT) and magnetic resonance 
imaging (MRI), for example. Three-dimensional (3D) images can be formed by 
stacking and interpolating between two-dimensional pictures produced from the 
scanning machines. Imaging an organ and visualizing its volume in three- 
dimensional space would be beneficial due to the lack of physical intrusion and the 
ease of data manipulation. However, the exploration of the three-dimensional 
volume image must be properly performed in order to fully exploit the advantages of 
virtually viewing an organ from the inside. 

With the progress of multi-detector computerized tomography (MDCT) and 
increasing temporal and spatial resolution of data sets, clinical use of computerized 
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tomographic angiography (CTA) is increasing. Vessel segmentation can be quite 
challenging, but it is needed to isolate vascular structures. 

Radiologists and other specialists have historically been trained to analyze 
scan data consisting of two-dimensional slices. However, while stacks of such slices 
5 may be useful for analysis, they do not provide an efficient or intuitive means to 
navigate through virtual organs, especially ones as complex as vascular structures. 
There remains a need for a virtual examination system providing data in a 
conventional format for analysis while, in addition, allowing an operator to easily 
navigate among vessels and vascular structures. 
10 SUMMARY 

These and other drawbacks and disadvantages of the prior art are addressed 
by a system and method for vessel segmentation using vesselness and edgeness. 

A system and corresponding method for vessel segmentation are provided. A 
system embodiment has an input adapter for receiving image data, a processor in 
15 signal communication with the input adapter, a pre-processing unit in signal 

communication with the processor for pre-processing the received image data, and a 
vessel segmentation unit in signal communication with the processor for segmenting 
vessels using pre-processed data. A corresponding method embodiment includes 
receiving image data, pre-processing the received data, and segmenting vessels 

i 

20 responsive to the pre-processed data. 
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These and other aspects, features and advantages of the present disclosure 
will become apparent from the following description of exemplary embodiments, 
which is to be read in connection with the accompanying drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 
5 Figure 1 shows a schematic diagram of a system for vessel segmentation 

using vesselness or edgeness in accordance with an embodiment of the present 
disclosure; 

Figure 2 shows a schematic diagram of a vessel model with a Gaussian 
luminance profile in accordance with an embodiment of the present disclosure; 
10 Figure 3 shows histograms from an aorta computerized tomographic 

angiography (CTA) scanning with and without contrast agent in accordance with an 
embodiment of the present disclosure; 

Figure 4 shows plot diagrams for CTA intensity ranges and the significant 
point in each range in accordance with an embodiment of the present disclosure; 
15 and 

Figure 5 shows a schematic diagram of a multi-level filter in accordance with 
an embodiment of the present disclosure. 

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 

The present disclosure provides system and method embodiments for vessel 
:0 segmentation using vesselness and edgeness. With the progress of multi-detector 
computerized tomography (MDCT) and increasing temporal and spatial resolution of 
data sets, clinical use of computerized tomographic angiography (CTA) data sets is 
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increasing. Vessel segmentation can be quite challenging, but it is needed to isolate 
vascular structures. 

As shown in Figure 1 , a system for vessel segmentation in CTA data sets 
using vesselness and edgeness, according to an illustrative embodiment of the 
present disclosure, is indicated generally by the reference numeral 100. The system 
100 includes at least one processor or Central processing unit (CPU) 102 in signal 
communication with a system bus 104. A read only memory (ROM) 106, a random 
access memory (RAM) 108, a display adapter 1 10, an I/O adapter 112, a user 
interface adapter 114, a communications adapter 128, and an imaging adapter 130 
are also in signal communication with the system bus 104. A display unit 1 16 is in 
signal communication with the system bus 104 via the display adapter 1 10. A disk 
storage unit 118, such as, for example, a magnetic or optical disk storage unit is in 
signal communication with the system bus 104 via the I/O adapter 1 12. A mouse 
120, and a keyboard 122 are in signal communication with the system bus 104 via 
the user interface adapter 114. An imaging device 132 is in signal communication 
with the system bus 104 via the imaging adapter 130. 

A pre-processing unit 170 and a vessel segmentation unit 180 are also 
included in the system 100 and in signal communication with the CPU 102 and the 
system bus 104. The exemplary pre-processing unit 170 includes a CTA pre-filtering 
portion, a multi-level vesselness computation portion, and an edgeness filter portion. 
The exemplary vessel segmentation unit 180 includes an integration portion for 
integrating vesselness and edgeness information. While the pre-processing unit 170 
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and the vessel segmentation unit 180 are illustrated as coupled to the at least one 
processor or CPU 102, these components are preferably embodied in computer 
program code stored in at least one of the memories 106, 108 and 118, wherein the 
computer program code is executed by the CPU 102. As will be recognized by 
those of ordinary skill in the pertinent art based on the teachings herein, alternate 
embodiments are possible, such as, for example, embodying some or all of the 
computer program code in registers located on the processor chip 102. Given the 
teachings of the disclosure provided herein, those of ordinary skill in the pertinent art 
will contemplate various alternate configurations and implementations of the pre- 
processing unit 170 and the multi-level vesselness computation unit 180, as well as 
the other elements of the system 100, while practicing within the scope and spirit of 
the present disclosure. 

i 

Preferred embodiments of the present disclosure provide a fast and robust 
method for CTA vessel segmentation using a two-step procedure, a pre-processing 
step and an interactive segmentation step. During the pre-processing step, a 
vesselness volume is computed by using a Hessian filter preceded by a CTA pre- 
filter. In order to enhance different size vessels, a Multum In Parvo (MIP) volume 
pyramid is created and multi-level vesselness is computed and merged. In addition, 
an edgeness volume, including the boundary information, is calculated and 
integrated with vesselness into an intermediate volume for vessel segmentation. At 
the interactive stage, a vessel central axis (VCA) is tracked and shown in real time by 
each click. After being initialized by VCAs, the system may take up to several 
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seconds to finish vessel segmentation. This method has been successfully 
implemented in CTA vessel extraction and evaluation due to its ease of use and 
reproducibility. 

CTA data sets are the exemplary context of the presently disclosed vessel 
5 segmentation method. The provided histogram-analysis based CTA pre-filter and 
fast vesselness computation method are adaptable to CTA datasets. Based on 
vesselness and edgeness, a fast front propagation method is presented to segment 
vessels with minimal user interaction. 

Embodiments of this disclosure focus on the following aspects of the presently 
10 disclosed method. A fast multi-level vesselness computation method is used. A 
preferred multi-level vesselness computation method uses a MlP-volume pyramid. 
In addition, it achieves the maximum response by filtering CTA data sets, called 
"CTA pre-filtering", before applying Hessian filtering. 

The Computation of Multilevel Vesselness is now described. A vessel is 
15 considered as a linear or tubular structure. Using a multi-scale line filter to detect or 
enhance tubular structures is one of the most popular methods, especially in 3D. 
Vesselness is a grade to assess tubular structures. The higher vesselness a voxel 
has, the greater the probability that it belongs to a tubular structure. Although not all 
of the tubular structures in vascular imaging are vessels, one can discriminate 
20 vessels from their surroundings with vesselness. 

Turning to Figure 2, a vessel model with a Gaussian luminance profile is 
indicated generally by the reference numeral 200. The vessel model 200 supposes 
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that a vessel model is a cylinder in which the Z-axis is the vessel central axis and the 
x-y plane is the vessel cross-section with Gaussian distribution, l 0 . 
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The Hessian matrix H is given by 
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The three eigenvalues and eigenvectors of the Hessian matrix H are 

i 

4 =0,1?! = (0,0,1); Eqn. 2 

K = ? 2 = (x 9 y,0) ; A, = v 3 - (-y,x,0) ; 

in which V 1} V 2 and V 3 are vectors perpendicular to each other. 



Based on this model, vesselness is scaled by the linearity of \ with /L, and A 39 
i.e., \ « 0 and ^ < \ « 0 , and the symmetry of \ and >\ , i.e., A,/^ ■« 1 . 

Different functions may be designated for S[i ne and S symm . For example, an 
15 exponent function or a polynomial function may be used. One issue with vesselness 
computations is that of the computational cost. The multi-scale convolutions of 2 nd 
derivatives can be very computational costly. The desire to reduce the 
computational cost without loss of vesselness quality is one of the main challenges in 
vesselness computation. 
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Experiments have been reported with methods using magnetic resonance 
angiography (MRA) data sets, including a head MRA case and a Carotid MRA case. 
In a later experiment, a bronchial case with threshold -180HU and a liver CTA case 
with threshold range of OHU ~ 300HU were discussed. Both cases had no adjacent 
bony structures and ignored calcium. Since some bones have tubular shapes, such 
as ribs compared to aorta, how to reduce the false-positive vesselness is another 
challenge to computing vesselness in CTA. 

Turning now to Figure 3, CTA Pre-filtering is now described with respect to 
histograms from an aorta CTA scanning with and without contrast agent, which are 
indicated generally by the reference numerals 300 and 350, respectively. The 
histograms 300 and 350 provide a comparison of an aorta CTA scanning with and 
without contrast agent, where the dashed line is the histogram scanned before 
injecting contrast agent, and the solid line the histogram scanned after injection. 

In CTA scanning, an iodinated contrast bolus is injected into a large vein and 
rapidly circulates through the heart and then reaches the arterial vessels. The 
histograms of Figure 3 show the effect of an IV contrast bolus in an abdominal CTA 
case. The histogram 300 has a range from -100HU to 900HU. The dashed line is 
the histogram before injecting contrast agent. The solid line is the result of contrast 
enhanced scanning. From the overview of both histograms, one can observe the 
shift of voxels in a low intensity range [-100, 0] to a high intensity range [50,500]. 
Thus, the range of Houndsfield units is markedly changed or enhanced. 
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An iodine contrast agent does help to distinguish vessels, but a vessel's 
intensity range still overlaps with other structures, especially lower-density bone and 
marrow. For example, when the vertebral artery goes through the cervical spines, or 
the anterior tibial artery is in close proximity to the tibia, they have a similar intensity 
5 range as bone surfaces such as low-density bone. Due to timing control and 
circulation of the contrast agent caused by plaques in vessels, the intensity of 
enhanced vessels is quite variable, from 100HU to 600HU. Because of this, a simple 
thresholding will not work well or consistently in CTA cases. 

The other obstacle of CTA vesselness computation is that in heavily diseased 

i 

10 arteries, calcium and other hard plaque adheres to the vessel wall and changes the 
local intensity profile, which makes it very difficult to compute the right eigenvalues 
and eigenvectors of the local Hessian matrix. Because high-density bone and 
calcification share the same intensity range, vesselness becomes discontinued or 
broken as a result of thresholding. Ideally, one would like to keep the calcium within 

15 vessels as parts of vessels segmented. However, a Hessian matrix will have 
difficulties in getting the right response when encountering hard plaques. 

An ideal vessel model is a Gaussian luminance profile, such as described by 
Equation 1 , where the highest intensity is at the middle of the tube and the intensity 
declines based on the Gaussian function at the boundary. However, this does not 

20 occur in CTA clinical practice, and one cannot directly apply the Hessian filter. 

Therefore, CTA data sets need a pre-filtering step to enhance the potential tubular 
structures before calculating vesselness. This CTA pre-filter should satisfy the 
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following requirements: Keep the Gaussian shape vessel luminal profile; adjust the 
volume intensity so that the maximum intensity within the vessels lumen becomes 
the maximum intensity of the volume; and normalize the intensity in order to compare 
the vesselness from different locations. 

As shown in Figure 4, three CTA intensity ranges and the Significant Point in 
each range are indicated generally by the reference numeral 400. A quadratic curve, 
called the Normalizing Roof Curve (NRC), is approximated by these three feature 
points. A look-up table based on NRC is calculated and used to map the original 
volume to keep the In-range voxels and dehance the Ex-vessel Low (ExL) range and 
Ex-vessel High (ExH) range voxels. 

Thus, the CTA histogram can be categorized into three ranges: Ex-vessel Low 
(ExL) range including air, fat and soft tissue; In-vessel (In) range including contrast 
enhanced vessel, low intensity bone and marrow etc.; and Ex-vessel High (ExH) 
range including bone and calcium. At both ends of each range, there is a Significant 
Point (SP). At this point, the slope changes significantly because the Region In 
Range (RIR) reduces considerably. SP is regarded as the statistical threshold to sort 
the different materials in the histogram. The SP is calculated by the intersection of 
two asymptotes, which are approached by tangent lines. In each range, one can 
calculate a pair of SP points A(SP 0 , SP^. From SP 0 the corresponding RIR starts 
growing massively. After SPi, the growing of RIR vanishes. 
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Since the histogram is separated into three ranges, three SP point pairs are 
calculated, i.e. Z\ E xl(SP 0 , SP,), A ]n (SP 0 , SP1J, 4 ExH (SPo, SP^.The pre-filter is set up 
as a roof-shape curve, specifically: 

Set SPi in An at the Peak Point 
5 Set SPq in A\ n at Left Verge Point 

Set SP 0 in ^ Ex h at Right Verge Point 

Multi-level filtering is now described. The linear filter is basically a multi-scale 
filter. That is, it filters with different radii (a) convolved at each voxel in the volume, 
and the maximum response is chosen as the vesselness at the current voxel. 

10 However, generally speaking, a vessel's radius can differ by as few as several voxels 
for a coronary artery, and up to a hundred voxels for an abdominal aorta. Therefore, 
a multi-scale filter must calculate all sizes in order to find the maximum response. 
On the other hand, large-scale filtering is very time consuming. 

Turning to Figure 5, a multi-level filter is indicated generally by the reference 

15 numeral 500. In a multi-level filter, a volume pyramid is created based on the original 
volume, called the Multum In Parvo (MIP) Volume. An MIP structure is widely used 
in computer graphics for 2D texture mapping, (MIP map). In an MIP Volume, Level 0 
is the original volume. The next level volume stores voxels recursively with half of 
the pre-level volume size, e.g., filtering or averaging over every 2x2x2 voxel. This 

20 process can reduce the volume size to one voxel or a certain level. With this volume 
pyramid, one can reconstruct any volume for which the size falls between two 
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integral levels. Every voxel can be interpolated with 16 nearest voxels to its neighbor 
integer level volumes, for 8 voxels on each level. 

Instead of using a multi-scale filter, a filter with fixed scale is used in multi- 
level filtering. The size of the filter is set up based on the smallest vessel expected in 
the original volume, such as, for example, a 5-voxel-diameter vessel. Then, the 
same filter is applied to compute the vesselness in other levels, such as 1.5, 2.0, 
etc. The resulting volume is scaled down to the 0-level size. Finally, results from all 
levels are merged into the 0-level volume by a maximum operator. 

V = max(V t ), 7 = 0,0.5,1.0,1.5, Eqn 3 

Calculating the vesselness using a multi-level filter avoids the time-consuming 
convolution of large-scale filters. The disadvantage is that it may introduce blurring 
in the resulting vesselness volume. Blurriness within vessels will not affect 
propagation or segmentation. If one can keep a sharp vesselness near a boundary, 
segmentation will get the same results as a multi-scale filter. In order to reduce the 
blurring near a boundary a gradient weight function is used before high-level 
vesselness is scaled down to 0-level. 

v /° = (1 -0 - gradient) * v, Eq n 4 

Front Propagation based on Vesselness is now described. Front propagation 
is an efficient fast marching level set method to track the monotonic advance of 
interfaces with a speed that does not change its sign, either expanding or shrinking. 
A front propagation is introduced in which the speed function is defined on 
vesselness. 
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Briefly, one can formalize the front propagation as the evolving function ip, 



namely: 



^+F|V^| = 0 Eqn. 5 

where F is the speed function. 

Letting Tbe defined as the implicit function such that ft={^+FV^ =0}, it is 



shown: 



d¥L = -FVy/ , where — = -FN 9 N = ^- x Eqn. 6 

dt dt 



Thus, the evolution of zero level-set of qj t equals the evolution of T t , a time- 
dependent implicit surface representing the evolving segmentation interface. 
10 Considering the special case of a monotonically advancing surface, such as 

for a front propagation, the speed function F is always positive and the normal vector 
N is always pointing outwards. Instead of numerical approximating the derivatives in 
Equation 5, one may explicit construct and update the front interface r t with Equation 
6. 

15 Timer-Tag Narrow Band is now described. The point to explicit construction of 

H is to create an active working zone, called a "narrow band", a local region around 
the front. The narrow band is constructed by active seeds. Each seed is defined by 
a vector of (P, t), of which P is the target position and t is the Timer. It explains that 
the current seed will take Mime to arrive at target position P. When t is positive, this 

20 seed is active. When t is less than or equal to zero, this seed is deactivated and 
merged into f t . 
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The narrow band is maintained by a heap data structure sorted by a timer f, 
called a "front-seed heap". The seed with smallest t is at the top of the heap. At 
each time step, the heap is checked and the front interface is marching outwards as 
below: 

For all seeds, decrease t 

For all the seeds of which t < 0: 

remove it from heap 

merge it into r t 

insert its neighbors into the heap 
Check the left seeds in heap of which t > 0: 

If P is in r t, remove the seed from the heap 

For each neighbor point NB of front point P, such as 6-neighbors, 18- 
neighbors, or 26-neighbors, the method checks only the outside A/6, i.e., 
outside the T, and initializes its timer with Equation 7: 

I LI 2 

t = j _[ — , where L = NB-P, F is the speed at point P Eqn. 7 

(Z ■ N)F 

The NB position and t are stored in the seed, and the seed is inserted into the 
front-seed heap. 

A speed function is now described, which defines the motion of the front. In 
order to segment vessels, a proper speed function is needed, which increases 
closing to the central part of vessels and vanishes at the boundary of vessels. Most 
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speed functions are defined directly on original images. A new speed function is now 
introduced that is defined based on vesselness. 

Two motions affect the speed function F in Equation 7, motion by curvature 
term and motion by the image term F/. 



.F = (1.0- £Ky)F I9 Eqn. 8 

where Ky is mean curvature andiT = V ( — — ) , e is the front smoothness 

| Vy/ 1 

control factor. 

Mostly, the F\ term is directly calculated from an original image by gradients, 
10 called static speed. Considering the vesselness, gradient and zero-crossing, a 
speed volume is initialized by using these two volumes. 

1.0 

~ ~h r. , 771 ~T\ (J- vesselness ^ ^ cons) Eqn. 9 

1.0 + C I grad 



where / con s is a constant term to keep the front moving. Near the zero- 
crossing points, this term is vanished. 
15 Since gradient and zero-crossing points are calculated while computing 

vesselness, instead of calculating F\ on fly, a speed volume is integrated and saved 

i 

directly after vesselness is calculated. By using the gradient and zero-crossing, the 
vesselness that might be blurred by the multi-level filter is converted to a sharpening 
speed at object boundaries. 
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Front Initialization is now described. The front is initialized by the vessel 
central axis (VGA) tracked in speed volume. In order to minimize user interaction, 
VCA is tracked based on the primary direction Vo in a Hessian matrix (see Equation 
2) from a user click. The VCAs along both primary directions (±Vo) are tracked. 
Unlike most ridge-tracking algorithms, preferred embodiments use a single-scale 
Hessian filter to estimate the primary direction, the same scale as used to calculate 
the vesselness. Advantages include: 

The vesselness or speed volume is already centralized, and single-scale is 
usually enough to track a VCA. 

VCA is used as an initial front to segment a vessel. Basically, if the initial front 
is located within a vessel, the results segmented by front propagation stay the same. 
That is to say, initial position is not critical to segmentation results. 

In addition, the VCA can be trimmed off when two central axes are close 
enough in space, that is, they can be merged into one central axis. The VCA always 
follows the directions of user clicks. At bifurcations, the user can control the vessel 
segmentation by selecting the interested branches. Since the VCA is tracked and 
displayed in real time, it gives the user an overview to interrogate the vessel that is 
going to be segmented. 

The exemplary method embodiment is integrated into an exemplary vascular 
measurement system and has been validated by CTA scannings from different parts 
of body and from different clinical institutions. 
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While exemplary embodiments have been described for use with CTA data 
sets, it will be recognized by those of ordinary skill in the pertinent art that alternate 
embodiments may be used for MRA, X-Ray Angiography (XRA) and Digital 
Subtraction Angiography (DSA) data sets, for example, even though MRA data sets 
may not show calcium deposits and bone. In addition, while the exemplary 
embodiments have been described as pre-processing values prior to segmentation, 
the values may alternately be computed during segmentation. In other words, the 
pre-processing and segmentation processes may be performed at the same time, 
wherein the vessel segmentation process directs the pre-processing of certain 
regions of the image volume data, as the user selects certain regions for 
segmentation (e.g., clicking on a region of a displayed image to segment a desired 
vessel or portion of the vessel). An advantage of this alternate embodiment is that 
only the values actually required for a selected segmentation would need to be 
computed for that segmentation. 

In other exemplary embodiments of the invention, a vessel segmentation and 
visualization system according to the invention enables selection and storage of 
multiple blood vessels for rapid reviewing at a subsequent time. For instance, a 
plurality of blood vessels that have been previously segmented, processed, 

„ i 

j 

annotated, etc. can be stored and later reviewed by selecting them one after another 
for rapid review. By way of further example, a plurality of different views may be 
simultaneously displayed in different windows (e.g., curved MPR, endoluminal view, 
etc.) for reviewing a selected blood vessel. When a user selects another stored (and 
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previously processed) blood vessel, all of the different views can be updated to 
include an image and relevant information associated with the newly selected blood 
vessel. In this manner, a user can select one or more multiple views that the user 
typically uses for reviewing blood vessels, for instance, and then selectively scroll 
5 through some or all of the stored blood vessels to have each of the views instantly 
updated with the selected blood vessels to rapidly review such stored set of vessels. 

The foregoing merely illustrates the principles of the disclosure. It will thus be 
appreciated that those skilled in the art will be able to devise numerous systems, 
apparatus and methods which, although not explicitly shown or described herein, 

10 embody the principles of the disclosure and are thus within the spirit and scope of 
the disclosure as defined by its Claims. 

For example, the methods and systems described herein could be applied to 
virtually examine an animal, fish or inanimate object. Besides the stated uses in the 
medical field, applications of the technique could be used to detect the contents of 

15 sealed objects that cannot be opened. The technique could also be used inside an 
architectural structure such as a building or cavern and enable the operator to 
navigate through the structure. 

These and other features and advantages of the present disclosure may be 
readily ascertained by one of ordinary skill in the pertinent art based on the teachings 

20 herein. It is to be understood that the teachings of the present disclosure may be 
implemented in various forms of hardware, software, firmware, special purpose 
processors, or combinations thereof. 
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Most preferably, the teachings of the present disclosure are implemented as a 
combination of hardware and software. Moreover, the software is preferably 
implemented as an application program tangibly embodied on a program storage 
unit. The application program may be uploaded to, and executed by, a machine 
comprising any suitable architecture. Preferably, the machine is implemented on a 
computer platform having hardware such as one or more central processing units 
(CPU), a random access memory (RAM), and input/output (I/O) interfaces. The 
computer platform may also include an operating system and microinstruction code. 
The various processes and functions described herein may be either part of the 
microinstruction code or part of the application program, or any combination thereof, 
which may be executed by a CPU. In addition, various other peripheral units may be 
connected to the computer platform such as an additional data storage unit and a 
printing unit. 

It is to be further understood that, because some of the constituent system 
components and methods depicted in the accompanying drawings are preferably 
implemented in software, the actual connections between the system components or 
the process function blocks may differ depending upon the manner in which 
embodiments of the present disclosure are programmed. Given the teachings 
herein, one of ordinary skill in the pertinent art will be able to contemplate these and 
similar implementations or configurations of the present invention. 

Although the illustrative embodiments have been described herein with 
reference to the accompanying drawings, it is to be understood that the present 
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invention is not limited to those precise embodiments, and that various changes and 
modifications may be effected therein by one of ordinary skill in the pertinent art 
without departing from the scope or spirit of the present disclosure. All such changes 
and modifications are intended to be included within the scope of the present 
5 invention as set forth in the appended Claims. 
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